!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine resistivity (bx_resis, by_resis, bz_resis, en_resis, ro, mx, my, mz, bx, by, bz, en, x, y, z, dx, dy, dz, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz
   real(8), intent(in) :: x(nx), y(ny), z(nz), dx, dy, dz
   real(8), dimension(nx, ny, nz), intent(inout) :: bx_resis, by_resis, bz_resis, en_resis, ro, mx, my, mz, bx, by, bz, en

   integer(4) :: ip1, jp1, kp1, im1, jm1, km1

   real(8) :: dx2, dy2, dz2, dxdx, dydy, dzdz
   real(8), dimension(nx, ny, nz) :: jx, jy, jz, fxen_resis, fyen_resis, fzen_resis, et

!----------------------------------------------------------------------------------------------------------------------------------|
!  initialization
!----------------------------------------------------------------------------------------------------------------------------------|
   dxdx = dx * dx
   dydy = dy * dy
   dzdz = dz * dz
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dz2 = dz * 2.0d0

   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      et(i, j, k) = 0.0d0
      jx(i, j, k) = 0.0d0
      jy(i, j, k) = 0.0d0
      jz(i, j, k) = 0.0d0
      bx_resis(i, j, k) = 0.0d0
      by_resis(i, j, k) = 0.0d0
      bz_resis(i, j, k) = 0.0d0
      en_resis(i, j, k) = 0.0d0
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  current density
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 2, nz - 1
      kp1 = k + 1
      km1 = k - 1
      do j = 2, ny - 1
         jp1 = j + 1
         jm1 = j - 1
         do i = 2, nx - 1
            ip1 = i + 1
            im1 = i - 1
            jx(i, j, k) = (bz(i, jp1, k) - bz(i, jm1, k)) / dy2 - (by(i, j, kp1) - by(i, j, km1)) / dz2
            jy(i, j, k) = (bx(i, j, kp1) - bx(i, j, km1)) / dz2 - (bz(ip1, j, k) - bz(im1, j, k)) / dx2
            jz(i, j, k) = (by(ip1, j, k) - by(im1, j, k)) / dx2 - (bx(i, jp1, k) - bx(i, jm1, k)) / dy2
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  resistivity
!----------------------------------------------------------------------------------------------------------------------------------|
   call model_resis(et, ro, mx, my, mz, bx, by, bz, en, jx, jy, jz, x, y, z, nx, ny, nz)

   max_et = max(maxval(et), max_et)

!----------------------------------------------------------------------------------------------------------------------------------|
!  variables update
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      fxen_resis(i, j, k) = et(i, j, k) * (by(i, j, k) * jz(i, j, k) - bz(i, j, k) * jy(i, j, k))
      fyen_resis(i, j, k) = et(i, j, k) * (bz(i, j, k) * jx(i, j, k) - bx(i, j, k) * jz(i, j, k))
      fzen_resis(i, j, k) = et(i, j, k) * (bx(i, j, k) * jy(i, j, k) - by(i, j, k) * jx(i, j, k))
   enddo
   enddo
   enddo

   do k = 2, nz - 1
      kp1 = k + 1
      km1 = k - 1
      do j = 2, ny - 1
         jp1 = j + 1
         jm1 = j - 1
         do i = 2, nx - 1
            ip1 = i + 1
            im1 = i - 1
            bx_resis(i, j, k) = (et(i, j, kp1) * jy(i, j, kp1) - et(i, j, km1) * jy(i, j, km1)) / dz2 -                            &
                                (et(i, jp1, k) * jz(i, jp1, k) - et(i, jm1, k) * jz(i, jm1, k)) / dy2
            by_resis(i, j, k) = (et(ip1, j, k) * jz(ip1, j, k) - et(im1, j, k) * jz(im1, j, k)) / dx2 -                            &
                                (et(i, j, kp1) * jx(i, j, kp1) - et(i, j, km1) * jx(i, j, km1)) / dz2
            bz_resis(i, j, k) = (et(i, jp1, k) * jx(i, jp1, k) - et(i, jm1, k) * jx(i, jm1, k)) / dy2 -                            &
                                (et(ip1, j, k) * jy(ip1, j, k) - et(im1, j, k) * jy(im1, j, k)) / dx2
            en_resis(i, j, k) = (fxen_resis(ip1, j, k) - fxen_resis(im1, j, k)) / dx2 +                                            &
                                (fyen_resis(i, jp1, k) - fyen_resis(i, jm1, k)) / dy2 +                                            &
                                (fzen_resis(i, j, kp1) - fzen_resis(i, j, km1)) / dz2
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine resistivity
